//kirk.cpp
#include <math.h>
//#include <stdio.h>
double besj1(double arg);
static double max,f1,f2;
//----------------------------------------------------------------------------
double mxf(double a,double b){
      double res=a;
      if(b>a) res=b;
return res;}
//----------------------------------------------------------------------------
void kirk0_init(double min1,double max1,double f10,double f20){
//***************** Parameter initialization **********************************
      f1=f10;
      f2=f20;
      max=mxf(fabs(min1),fabs(max1));
return;}
//-----------------------------------------------------------------------------
double kirk0(double u,double v,double a){
//********* Phase function computation at single point ************************
      double bessel0(double);
      double psi,psi0=0,pi=3.1415926;
      if(a<0) psi0=pi;
      psi=psi0+bessel0(fabs(a)/max)*cos(2*pi*(f1*u+f2*v));
return psi;}
void kirk1_init(double min1,double max1,double f10,double f20){
//***************** Parameter initialization **********************************
      f1=f10;
      f2=f20;
      max=mxf(fabs(min1),fabs(max1));
return;}
//-----------------------------------------------------------------------------
double kirk1(double u,double v,double a){
//********* Phase function computation at single point ************************
      double bessel1(double);
      double psi,psi0=0,pi=3.1415926;
      if(a<0) psi0=pi;
      psi=psi0+bessel1(besj1(2*M_PI)*fabs(a)/max)*cos(2*pi*(f1*u+f2*v));
return psi;}
double bessel0( double arg) {
//************ Transcendent bessel equation solution **************************
double arg1;
arg1=arg;
if(arg>1) arg1=1;if(arg<0) arg1=0;
double a[101]=
{2.38, 2.37, 2.355, 2.335, 2.32, 2.3, 2.285, 2.265, 2.25, 2.23, 2.21, 2.19, 2.175, 2.16,
2.145, 2.13, 2.105, 2.085, 2.075, 2.05, 2.03, 2.02, 2.015, 1.985, 1.97, 1.955, 1.935,
1.915, 1.91, 1.885, 1.86, 1.85, 1.835, 1.81, 1.795, 1.78, 1.765, 1.755, 1.73, 1.715, 1.69,
1.675, 1.66, 1.645, 1.625, 1.61, 1.59, 1.57, 1.56, 1.54, 1.52, 1.5, 1.48, 1.46, 1.45, 1.43,
1.41, 1.39, 1.37, 1.35, 1.33, 1.32, 1.3, 1.28, 1.26, 1.24, 1.22, 1.2, 1.18,
1.16, 1.14, 1.12, 1.09, 1.07, 1.05, 1.03, 1.01, 0.985, 0.965, 0.94,
0.91, 0.89, 0.86, 0.84, 0.8, 0.79, 0.76, 0.73, 0.7, 0.67, 0.64, 0.605,
0.57, 0.535, 0.49, 0.45, 0.4, 0.345, 0.285, 0.2,0.000};
int index;
index=(int)(arg1*100);
return(a[index]);}
double besj1(double arg);
double bessel1(double arg){
//*** The solution of the Bessel-1 transcendent equation by dihotomia ****
float a,b,c;
if(0>arg)  return  0;
if(M_PI_2<arg) return besj1(M_PI_2);
a=0;
b=1.84;
nach: c=(a+b)/2;
      if(fabs(besj1(c)-arg)<0.02) goto con;
      if((besj1(c)-arg)<0) a=c; else b=c;
      goto nach;
con:  return c;}
